---
output: 
  pdf_document:
    citation_package: natbib
    keep_tex: true
    fig_caption: true
    latex_engine: pdflatex
    template: ~/Dropbox/Software/templates/svm-r-markdown-templates-master/svm-latex-ms.tex
# for template details, see http://svmiller.com/blog/2016/02/svm-r-markdown-manuscript/    
title: "Appendix to Decomposing Public Opinion Variation"
# thanks: "**Current version**: `r format(Sys.time(), '%B %d, %Y')`; **Corresponding author**: b.e.lauderdale@lse.ac.uk"  
# author:
# - name: Benjamin E Lauderdale
#   affiliation: London School of Economics and Political Science
# - name: Chris Hanretty
#   affiliation: Royal Holloway, University of London
# - name: Nick Vivyan
#   affiliation: Durham University
# abstract: "In this appendix, we report decompositions from a series of robustness checks mentioned in our article."
keywords: ""
date: "`r format(Sys.time(), '%B %d, %Y')`"
geometry: margin=1in
fontfamily: mathpazo
fontsize: 10pt
# spacing: double
# bibliography: master.bib
# biblio-style: apsr
---

```{r, include=FALSE}

# Libraries

library(knitr)
library(pander)
library(xtable)
library(rstan)
library(abind)
library(ggplot2)
library(ggtern)
library(ggrepel)
library(MCMCpack)
library(reshape2)
library(dplyr)


# Settings

set.seed(42)
REFIT <- FALSE # Set to TRUE to fit models, FALSE to generate report from saved models
HMCSample <- 2000
HMCBurn <- 500
HMCChains <- 2
rstan_options(auto_write = TRUE)
options(mc.cores = HMCChains)

# Function definitions

rmse <- function(x1,x2) sqrt(mean((x1-x2)^2))
extract_chain <- function(var) extract(posterior,var)[[1]]
fadecol <- function(colours,alpha=1){
  rgb.dat <- col2rgb(colours)
  return(rgb(rgb.dat[1,],rgb.dat[2,],rgb.dat[3,],alpha*255,maxColorValue=255))
}  

```


```{r, include=FALSE}

# load survey data
alldata <- read.csv("ssi-spatial-all-merged-with-wave1.tab",sep="\t")
nobs_full <- nrow(alldata)

# Select only respondents interviewed twice...

panelobs <- rowSums(!is.na(alldata[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")])) != 0

t.test(alldata$pkscale ~ panelobs)
alldata <- subset(alldata,subset=panelobs)
nobs_panel <- nrow(alldata)

alldata$pkscale2 <- (alldata$pkscale - min(alldata$pkscale))/(max(alldata$pkscale) - min(alldata$pkscale))
varnames <- c("Marijuana","Environment","Social_Security","Guns","Health","Immigration","Taxes","Abortion","Medicare","Gay_Rights","Unions","Contraception","Education")
alldata$partycol <- c("#0000FF","#0000DD","#0000BB","#BB00BB","#BB0000","#DD0000","#FF0000")[alldata$partyid]

```

```{r collapse, include = FALSE}

w1qs <- c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
          "eq_contraception","eq_education")
w2qs <- paste0(w1qs, "_POST")

alldata.coll <- alldata
### Collapse categories 1 and 2, 6 and 7
for (v in c(w1qs, w2qs)) {
    alldata.coll[,v] <- dplyr::recode(alldata.coll[,v],
                                 `1` = 1L,
                                 `2` = 1L,
                                 `3` = 2L,
                                 `4` = 3L,
                                 `5` = 4L,
                                 `6` = 5L,
                                 `7` = 5L)
}

# Stan Model Data
Y1 <- alldata.coll[,c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_contraception","eq_education")]
for (j in 1:ncol(Y1)) Y1[,j] <- as.integer(Y1[,j])
Y1 <- as.matrix(Y1)
colnames(Y1) <- varnames
Y1 <- replace(Y1,is.na(Y1),-1L)

Y2 <- alldata.coll[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")]
for (j in 1:ncol(Y2)) Y2[,j] <- as.integer(Y2[,j])
Y2 <- as.matrix(Y2)
colnames(Y2) <- varnames
Y2 <- replace(Y2,is.na(Y2),-1L)

Ycoll5 <- abind(Y1,Y2,along=3)

alldata.coll <- alldata
### Collapse 1 and 2; 3 and 4 and 5; 6 and 7
for (v in c(w1qs, w2qs)) {
    alldata.coll[,v] <- dplyr::recode(alldata.coll[,v],
                                 `1` = 1L,
                                 `2` = 1L,
                                 `3` = 2L,
                                 `4` = 2L,
                                 `5` = 2L,
                                 `6` = 3L,
                                 `7` = 3L)
}

# Stan Model Data
Y1 <- alldata.coll[,c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_contraception","eq_education")]
for (j in 1:ncol(Y1)) Y1[,j] <- as.integer(Y1[,j])
Y1 <- as.matrix(Y1)
colnames(Y1) <- varnames
Y1 <- replace(Y1,is.na(Y1),-1L)

Y2 <- alldata.coll[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")]
for (j in 1:ncol(Y2)) Y2[,j] <- as.integer(Y2[,j])
Y2 <- as.matrix(Y2)
colnames(Y2) <- varnames
Y2 <- replace(Y2,is.na(Y2),-1L)

Ycoll3 <- abind(Y1,Y2,along=3)

alldata.coll <- alldata
### Collapse 1, 2 and 3; 4; 5, 6 and 7
for (v in c(w1qs, w2qs)) {
    alldata.coll[,v] <- dplyr::recode(alldata.coll[,v],
                                 `1` = 1L,
                                 `2` = 1L,
                                 `3` = 1L,
                                 `4` = 2L,
                                 `5` = 3L,
                                 `6` = 3L,
                                 `7` = 3L)
}

# Stan Model Data
Y1 <- alldata.coll[,c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_contraception","eq_education")]
for (j in 1:ncol(Y1)) Y1[,j] <- as.integer(Y1[,j])
Y1 <- as.matrix(Y1)
colnames(Y1) <- varnames
Y1 <- replace(Y1,is.na(Y1),-1L)

Y2 <- alldata.coll[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")]
for (j in 1:ncol(Y2)) Y2[,j] <- as.integer(Y2[,j])
Y2 <- as.matrix(Y2)
colnames(Y2) <- varnames
Y2 <- replace(Y2,is.na(Y2),-1L)

Ycoll3b <- abind(Y1,Y2,along=3)

# Stan Model Data
Y1 <- alldata[,c("eq_marijuana","eq_environment","eq_socialsecurity",
                "eq_guns","eq_health","eq_immigration","eq_taxes",
                "eq_abortion","eq_medicare","eq_gays","eq_unions",
                "eq_contraception","eq_education")]
for (j in 1:ncol(Y1)) Y1[,j] <- as.integer(Y1[,j])
Y1 <- as.matrix(Y1)
colnames(Y1) <- varnames
Y1 <- replace(Y1,is.na(Y1),-1L)

Y2 <- alldata[,c("eq_marijuana_POST","eq_environment_POST","eq_socialsecurity_POST",
                 "eq_guns_POST","eq_health_POST","eq_immigration_POST","eq_taxes_POST",
                 "eq_abortion_POST","eq_medicare_POST","eq_gays_POST","eq_unions_POST",
                 "eq_contraception_POST","eq_education_POST")]
for (j in 1:ncol(Y2)) Y2[,j] <- as.integer(Y2[,j])
Y2 <- as.matrix(Y2)
colnames(Y2) <- varnames
Y2 <- replace(Y2,is.na(Y2),-1L)

Y <- abind(Y1,Y2,along=3)
```

```{r fakedata, include = FALSE}

load("stanout.Rdata")  
posterior <- posterior1
alpha_est <- apply(extract_chain("alpha"),c(2,3),mean)
beta_est <- apply(extract_chain("beta"),c(2),mean)
sigma_est <- apply(extract_chain("sigma"),c(2),mean)
omega_est <- apply(extract_chain("omega"),c(2),mean)

fakef <- function(Y,alpha,beta,sigma,omega) {

    nRespondents <- dim(Y)[1]
    nItems <- dim(Y)[2]
    nPeriods <- dim(Y)[3]
    
### Generate ideologies
    theta <- rnorm(n = nRespondents, mean = 0, sd = 1)
    
### Generate idiosyncracies per issue per respondent
    
    nu_std <- matrix(rnorm(n = nItems, mean = 0,
                   sd = 1),
                  nrow=nRespondents,
                  ncol=nItems)
    nu <- nu_std
    for (j in 1:nItems) nu[,j] <- nu_std[,j]*omega[j]
    
### Go through each cell and replace real data with fake data    
    
    Ysim <- Y
    for (i in 1:nRespondents){
        for (j in 1:nItems){
            for (p in 1:2){
                if (Ysim[i,j,p] != -1) Ysim[i,j,p] <- 1 + sum((beta[j]*theta[i] + nu[i,j] + rnorm(1)*sigma[j] - alpha[j,]) > 0)
            }}}
    
    return(Ysim)
}

fakeY1 <- fakef(Y,alpha_est,beta_est,sigma_est,omega_est)
fakeY2 <- fakef(Y,alpha_est,beta_est,sqrt(sigma_est^2+omega_est^2),rep(0,length(omega_est)))

```

```{r, include=FALSE}

standata <- list(
  Y = Y, # Responses to 5 pt scales
  N_respondents = nrow(Y),
  N_items = ncol(Y),
  N_categories = max(Y)
)

standata5 <- list(
  Y = Ycoll5, # Responses to 5 pt scales
  N_respondents = nrow(Ycoll5),
  N_items = ncol(Ycoll5),
  N_categories = max(Ycoll5)
)

standata3 <- list(
  Y = Ycoll3, # Responses to 3 pt scales
  N_respondents = nrow(Ycoll3),
  N_items = ncol(Ycoll3),
  N_categories = max(Ycoll3)
)

standata3b <- list(
  Y = Ycoll3b, # Responses to 3 pt scales
  N_respondents = nrow(Ycoll3b),
  N_items = ncol(Ycoll3b),
  N_categories = max(Ycoll3b)
)


fakedata1 <- list(
  Y = fakeY1, # Responses to 7 pt scales
  N_respondents = nrow(fakeY1),
  N_items = ncol(fakeY1),
  N_categories = max(fakeY1)
)

fakedata2 <- list(
  Y = fakeY2, # Responses to 7 pt scales
  N_respondents = nrow(fakeY2),
  N_items = ncol(fakeY2),
  N_categories = max(fakeY2)
)

## Model with spatial dimension

stancode1 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std;
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta;
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects

}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] = sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[j]*beta_std[j];
    sigma[j] = sigma_std[j]/sqrt(norm[j]); 
    omega[j] = omega_std[j]/sqrt(norm[j]); 
    beta[j] = beta_std[j]/sqrt(norm[j]);
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] = nu_std[i,j]*omega[j];
      eta[i,j] = beta[j]*theta[i] + nu[i,j]; // spatial model
      pi[i,j,1] = 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] = Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] = Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std ~ normal(0,1);
  theta ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] = beta[j]*beta[j];
    decomposition[j,2] = omega[j]*omega[j];
    decomposition[j,3] = sigma[j]*sigma[j];
  }

}

'

## Model without idiosyncracy

stancode2 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std;
  vector<lower=0>[N_items] sigma_std;
  vector[N_respondents] theta;
  ordered[N_categories-1] alpha[N_items];
  
}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real eta[N_respondents,N_items];
  vector[N_categories] pi[N_respondents,N_items];

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] = sigma_std[j]*sigma_std[j] + beta_std[j]*beta_std[j];
    sigma[j] = sigma_std[j]/sqrt(norm[j]); 
    omega[j] = 0.0;
    beta[j] = beta_std[j]/sqrt(norm[j]);
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      eta[i,j] = beta[j]*theta[i]; // spatial model
      pi[i,j,1] = 1.0 - Phi((eta[i,j] - alpha[j,1])/sigma[j]);
      for (k in 2:(N_categories-1)) pi[i,j,k] = Phi((eta[i,j] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j] - alpha[j,k])/sigma[j]);
      pi[i,j,N_categories] = Phi((eta[i,j] - alpha[j,N_categories-1])/sigma[j]);
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j]);
    }
  }
  sigma_std ~ normal(0,1);
  beta_std ~ normal(0,1);
  theta ~ normal(0,1);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] = beta[j]*beta[j];
    decomposition[j,2] = omega[j]*omega[j];
    decomposition[j,3] = sigma[j]*sigma[j];
  }

}

'

## Model with changing spatial position

stancode3 <- '

data {
  int<lower = 1> N_respondents; // Number of respondents 
  int<lower = 1> N_items; // Number of questions
  int<lower = 1> N_categories; // Number of response categories
  int<lower = -1,upper=N_categories> Y[N_respondents,N_items,2]; // Array of response data  
}

parameters {
  vector[N_items] beta_std;
  vector<lower=0>[N_items] omega_std;
  vector<lower=0>[N_items] sigma_std;
  vector[2] theta[N_respondents];
  ordered[N_categories-1] alpha[N_items];
  real nu_std[N_respondents,N_items]; // person-item random effects
  real<lower=-1,upper=1> rho;
}

transformed parameters {

  real<lower=0> norm[N_items];
  real<lower=0> omega[N_items];
  real<lower=0> sigma[N_items];
  vector[N_items] beta;
  real nu[N_respondents,N_items]; // person-item random effects
  real eta[N_respondents,N_items,2];
  vector[N_categories] pi[N_respondents,N_items,2];
  vector[2] theta_prior_mean;
  matrix[2,2] theta_prior_variance;

  theta_prior_mean[1] = 0;
  theta_prior_mean[2] = 0;
  theta_prior_variance[1,1] = 1;
  theta_prior_variance[1,2] = rho;
  theta_prior_variance[2,1] = rho;
  theta_prior_variance[2,2] = 1;

  for(j in 1:N_items) {  // Normalize so that latent scale has variance equal to one for each item
    norm[j] = sigma_std[j]*sigma_std[j] + omega_std[j]*omega_std[j] + beta_std[j]*beta_std[j];
    sigma[j] = sigma_std[j]/sqrt(norm[j]); 
    omega[j] = omega_std[j]/sqrt(norm[j]); 
    beta[j] = beta_std[j]/sqrt(norm[j]);
  }

  for(i in 1:N_respondents){
    for(j in 1:N_items){
      nu[i,j] = nu_std[i,j]*omega[j];
      for (p in 1:2){
        eta[i,j,p] = beta[j]*theta[i,p] + nu[i,j]; // spatial model
        pi[i,j,p,1] = 1.0 - Phi((eta[i,j,p] - alpha[j,1])/sigma[j]);
        for (k in 2:(N_categories-1)) pi[i,j,p,k] = Phi((eta[i,j,p] - alpha[j,k-1])/sigma[j]) - Phi((eta[i,j,p] - alpha[j,k])/sigma[j]);
        pi[i,j,p,N_categories] = Phi((eta[i,j,p] - alpha[j,N_categories-1])/sigma[j]);
      }
    }
  }
}

model {

  for(j in 1:N_items){
    for(i in 1:N_respondents){
      nu_std[i,j] ~ normal(0,1);
      if (Y[i,j,1] != -1) Y[i,j,1] ~ categorical(pi[i,j,1]);
      if (Y[i,j,2] != -1) Y[i,j,2] ~ categorical(pi[i,j,2]);
    }
  }
  omega_std ~ normal(0,1);
  sigma_std ~ normal(0,1);
  beta_std ~ normal(0,1);
  theta ~ multi_normal(theta_prior_mean,theta_prior_variance);
}

generated quantities {

  simplex[3] decomposition[N_items];
  for(j in 1:N_items) {
    decomposition[j,1] = beta[j]*beta[j];
    decomposition[j,2] = omega[j]*omega[j];
    decomposition[j,3] = sigma[j]*sigma[j];
  }

}

'

geninit3 <- function() return(list(
  beta_std = rep(0,standata3$N_items),
  sigma_std = rep(1,standata3$N_items),
  omega_std = rep(1,standata3$N_items),  
  theta = rnorm(standata3$N_respondents,0,1),
  nu_std = matrix(rnorm(standata3$N_respondents*standata3$N_items,0,0.1),standata3$N_respondents,standata3$N_items),
  alpha = matrix(seq(-2,2,length.out=standata3$N_categories-1),standata3$N_items,standata3$N_categories-1,byrow=TRUE)
))

geninit5 <- function() return(list(
  beta_std = rep(0,standata5$N_items),
  sigma_std = rep(1,standata5$N_items),
  omega_std = rep(1,standata5$N_items),  
  theta = rnorm(standata5$N_respondents,0,1),
  nu_std = matrix(rnorm(standata5$N_respondents*standata5$N_items,0,0.1),standata5$N_respondents,standata5$N_items),
  alpha = matrix(seq(-2,2,length.out=standata5$N_categories-1),standata5$N_items,standata5$N_categories-1,byrow=TRUE)
))

geninitF <- function() return(list(
  beta_std = rep(0,fakedata1$N_items),
  sigma_std = rep(1,fakedata1$N_items),
  omega_std = rep(1,fakedata1$N_items),  
  theta = rnorm(fakedata1$N_respondents,0,1),
  nu_std = matrix(rnorm(fakedata1$N_respondents*fakedata1$N_items,0,0.1),fakedata1$N_respondents,fakedata1$N_items),
  alpha = matrix(seq(-2,2,length.out=fakedata1$N_categories-1),fakedata1$N_items,fakedata1$N_categories-1,byrow=TRUE)
))


geninitI <- function() return(list(
  beta_std = rep(0,standata$N_items),
  sigma_std = rep(1,standata$N_items),
  ## omega_std = rep(1,standata$N_items),  
  theta = rnorm(standata$N_respondents,0,1),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

geninitC <- function() return(list(
  beta_std = rep(0,standata$N_items),
  sigma_std = rep(1,standata$N_items),
  omega_std = rep(1,standata$N_items),  
  theta = cbind(rnorm(standata$N_respondents,0,1),rnorm(standata$N_respondents,0,1)),
  nu_std = matrix(rnorm(standata$N_respondents*standata$N_items,0,0.1),standata$N_respondents,standata$N_items),
  alpha = matrix(seq(-2,2,length.out=standata$N_categories-1),standata$N_items,standata$N_categories-1,byrow=TRUE)
))

if (REFIT) {
 
    posterior3Category <- stan(
        model_code = stancode1, 
        data = standata3,
        init= geninit3,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )    
  
    posterior3CategoryB <- stan(
        model_code = stancode1, 
        data = standata3b,
        init= geninit3,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )    
          
    posterior5Category <- stan(
        model_code = stancode1, 
        data = standata5,
        init= geninit5,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )
    
    posteriorFakeData1 <- stan(
        model_code = stancode1, 
        data = fakedata1,
        init= geninitF,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )    
    
    posteriorFakeData2 <- stan(
        model_code = stancode1, 
        data = fakedata2,
        init= geninitF,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )
    
    posteriorNoIdiosyncrasy <- stan(
        model_code = stancode2, 
        data = standata,
        init= geninitI,
        pars=c(
            "alpha","beta","theta","sigma","omega","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )
    
    posteriorIdeologyChange <- stan(
        model_code = stancode3, 
        data = standata,
        init= geninitC,
        pars=c(
            "alpha","beta","theta","sigma","omega","rho","decomposition"
        ), 
        iter = HMCSample + HMCBurn,
        warmup = HMCBurn, 
        chains = HMCChains,
        refresh = 1,
        save_dso = FALSE
    )
    
        
  save(posterior3Category,posterior3CategoryB,posterior5Category, posteriorFakeData1,posteriorFakeData2, posteriorNoIdiosyncrasy,posteriorIdeologyChange, file = "stanout_appendix.Rdata")
  
} else load("stanout_appendix.Rdata")  

load("stanout.Rdata")  


```

```{r fakeanalysis, include = FALSE}

posterior <- posterior0
decomposition_chain_0 <- extract_chain("decomposition")
decomposition_est_0 <- t(apply(decomposition_chain_0,c(2,3),mean))
colnames(decomposition_est_0) <- varnames
rownames(decomposition_est_0) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0)

posterior <- posterior0bcit
decomposition_chain_0bcit <- extract_chain("decomposition")
decomposition_est_0bcit <- t(apply(decomposition_chain_0bcit,c(2,3),mean))
rownames(decomposition_est_0bcit) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0bcit)

posterior <- posterior0bleg
decomposition_chain_0bleg <- extract_chain("decomposition")
decomposition_est_0bleg <- t(apply(decomposition_chain_0bleg,c(2,3),mean))
rownames(decomposition_est_0bleg) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_0bleg)

posterior <- posterior1
decomposition_chain_1 <- extract_chain("decomposition")
decomposition_est_1 <- t(apply(decomposition_chain_1,c(2,3),mean))
decomposition_se_1 <- t(apply(decomposition_chain_1,c(2,3),sd))
colnames(decomposition_est_1) <- colnames(decomposition_se_1) <- varnames
rownames(decomposition_est_1) <- rownames(decomposition_se_1) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1)

posterior <- posterior1b
decomposition_chain_1b <- extract_chain("decomposition")
decomposition_est_1b <- t(apply(decomposition_chain_1b,c(2,3),mean))
decomposition_se_1b <- t(apply(decomposition_chain_1b,c(2,3),sd))
colnames(decomposition_est_1b) <- colnames(decomposition_se_1b) <- varnames
rownames(decomposition_est_1b) <- rownames(decomposition_se_1b) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1b)

posterior <- posterior1c
decomposition_chain_1c <- extract_chain("decomposition")
decomposition_est_1c <- t(apply(decomposition_chain_1c,c(2,3),mean))
decomposition_se_1c <- t(apply(decomposition_chain_1c,c(2,3),sd))
colnames(decomposition_est_1c) <- colnames(decomposition_se_1c) <- varnames
rownames(decomposition_est_1c) <- rownames(decomposition_se_1c) <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_1c)

posterior <- posterior2
decomposition_chain_2 <- extract_chain("decomposition")
decomposition_est_2 <- apply(decomposition_chain_2,c(2,3,4),mean)
decomposition_se_2 <- apply(decomposition_chain_2,c(2,3,4),sd)
dimnames(decomposition_est_2)[[1]] <- dimnames(decomposition_se_2)[[1]] <- varnames
dimnames(decomposition_est_2)[[2]] <- dimnames(decomposition_se_2)[[2]] <- c("Low Knowledge","High Knowledge")
dimnames(decomposition_est_2)[[3]] <- dimnames(decomposition_se_2)[[3]] <- c("Ideology","Idiosyncrasy","Instability")
print(decomposition_est_2)

### Get the decomposition for the fake data analysis
posterior <- posteriorFakeData1
decomp_chain_fake1 <- extract_chain("decomposition")
decomp_est_fake1 <- t(apply(decomp_chain_fake1,c(2,3),mean))
colnames(decomp_est_fake1) <- varnames
rownames(decomp_est_fake1) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_fake1)
rowMeans(decomp_est_fake1)

### Get the decomposition for the fake data analysis
posterior <- posteriorFakeData2
decomp_chain_fake2 <- extract_chain("decomposition")
decomp_est_fake2 <- t(apply(decomp_chain_fake2,c(2,3),mean))
colnames(decomp_est_fake2) <- varnames
rownames(decomp_est_fake2) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_fake2)
rowMeans(decomp_est_fake2)

### Get the decomp for the collapsed to 5 categories analysis
posterior <- posterior5Category
decomp_chain_cut5 <- extract_chain("decomposition")
decomp_est_cut5 <- t(apply(decomp_chain_cut5,c(2,3),mean))
colnames(decomp_est_cut5) <- varnames
rownames(decomp_est_cut5) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_cut5)
rowMeans(decomp_est_cut5)

### Get the decomp for the collapsed to 3 categories analysis
posterior <- posterior3Category
decomp_chain_cut3 <- extract_chain("decomposition")
decomp_est_cut3 <- t(apply(decomp_chain_cut3,c(2,3),mean))
colnames(decomp_est_cut3) <- varnames
rownames(decomp_est_cut3) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_cut3)
rowMeans(decomp_est_cut3)

### Get the decomp for the collapsed to 3 categories analysis, version B
posterior <- posterior3CategoryB
decomp_chain_cut3b <- extract_chain("decomposition")
decomp_est_cut3b <- t(apply(decomp_chain_cut3b,c(2,3),mean))
colnames(decomp_est_cut3b) <- varnames
rownames(decomp_est_cut3b) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_cut3b)
rowMeans(decomp_est_cut3b)

### Get the decomp for the ideology alone analysis
posterior <- posteriorNoIdiosyncrasy
decomp_chain_ideo <- extract_chain("decomposition")
decomp_est_ideo <- t(apply(decomp_chain_ideo,c(2,3),mean))
colnames(decomp_est_ideo) <- varnames
rownames(decomp_est_ideo) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_ideo)
rowMeans(decomp_est_ideo)

### Get the decomp for the ideology alone analysis
posterior <- posteriorIdeologyChange
decomp_chain_change <- extract_chain("decomposition")
decomp_est_change <- t(apply(decomp_chain_change,c(2,3),mean))
colnames(decomp_est_change) <- varnames
rownames(decomp_est_change) <- c("Ideology","Idiosyncrasy","Instability")
print(decomp_est_change)
rowMeans(decomp_est_change)
rho_chain <- extract_chain("rho")

decomp_avg_chain_cut5 <- apply(decomp_chain_cut5,c(1,3),mean)
decomp_avg_cut5_int <- 100*apply(decomp_avg_chain_cut5,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

decomp_avg_chain_cut3 <- apply(decomp_chain_cut3,c(1,3),mean)
decomp_avg_cut3_int <- 100*apply(decomp_avg_chain_cut3,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

decomp_avg_chain_cut3b <- apply(decomp_chain_cut3b,c(1,3),mean)
decomp_avg_cut3b_int <- 100*apply(decomp_avg_chain_cut3b,
                                           2,
                                           quantile,
                                           c(0.025,0.975))


decomp_avg_chain_fake1 <- apply(decomp_chain_fake1,c(1,3),mean)
decomp_avg_fake1_int <- 100*apply(decomp_avg_chain_fake1,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

decomp_avg_chain_fake2 <- apply(decomp_chain_fake2,c(1,3),mean)
decomp_avg_fake2_int <- 100*apply(decomp_avg_chain_fake2,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

decomp_avg_chain_ideo <- apply(decomp_chain_ideo,c(1,3),mean)
decomp_avg_ideo_int <- 100*apply(decomp_avg_chain_ideo,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

decomp_avg_chain_change <- apply(decomp_chain_change,c(1,3),mean)
decomp_avg_change_int <- 100*apply(decomp_avg_chain_change,
                                           2,
                                           quantile,
                                           c(0.025,0.975))

```

```{r, include = FALSE}

decomptable <- rbind(
apply(decomposition_chain_1,c(3),mean),
apply(decomposition_chain_0,c(3),mean),
apply(decomp_chain_ideo,c(3),mean),
apply(decomposition_chain_0bcit,c(3),mean),
apply(decomposition_chain_0bleg,c(3),mean),
apply(decomposition_chain_1b,c(3),mean),
apply(decomposition_chain_1c,c(3),mean),
apply(decomp_chain_change,c(3),mean),
apply(decomposition_chain_2[,,1,],c(3),mean),
apply(decomposition_chain_2[,,2,],c(3),mean),
apply(decomp_chain_fake1,c(3),mean),
apply(decomp_chain_fake2,c(3),mean),
apply(decomp_chain_cut5,c(3),mean),
apply(decomp_chain_cut3,c(3),mean),
apply(decomp_chain_cut3b,c(3),mean))

get_range <- function(chain, slice = 1) {
    if (length(dim(chain)) == 3) { 
        retval <- apply(chain, c(1, 3), mean)
        retval <- apply(retval, 2, quantile, c(0.025, 0.975))
        retval[1,] <- floor(100 * retval[1, ])
        retval[2,] <- ceiling(100 * retval[2, ])
        retval <- apply(retval, 2, function(x) paste0("(", x[1], " - ", x[2], ")", collapse = ""))
    } else {
        retval <- apply(chain, c(1, 3, 4), mean)
        retval <- retval[,slice,]
        retval <- apply(retval, 2, quantile, c(0.025, 0.975))
        retval[1,] <- floor(100 * retval[1, ])
        retval[2,] <- ceiling(100 * retval[2, ])
        retval <- apply(retval, 2, function(x) paste0("(", x[1], " - ", x[2], ")", collapse = ""))
    }
}


inttabl <- rbind(
    get_range(decomposition_chain_1),
    get_range(decomposition_chain_0),
    get_range(decomp_chain_ideo),
    get_range(decomposition_chain_0bcit),
    get_range(decomposition_chain_0bleg),
    get_range(decomposition_chain_1b),
    get_range(decomposition_chain_1c),
    get_range(decomp_chain_change),
    get_range(decomposition_chain_2, slice = 1),
    get_range(decomposition_chain_2, slice = 2),
    get_range(decomp_chain_fake1),
    get_range(decomp_chain_fake2),
    get_range(decomp_chain_cut5),
    get_range(decomp_chain_cut3),
    get_range(decomp_chain_cut3b))

variantnames <- c("Main model",
  "No ideology",
  "No idiosyncrasy",
  "One wave",
  "One wave, state legislator survey",
  "2 ideology dimensions",
  "3 ideology dimensions",
  "Ideology by period",
  "Low knowledge respondents",
  "High knowledge respondents",
  "Simulated data, estimated parameters",
  "Simulated data, no idiosyncrasy",
  "Responses grouped 12-3-4-5-67",
  "Responses grouped 12-345-67",
  "Responses grouped 123-4-567"
)

outtab <- matrix(data = "&nbsp;",
                 nrow = length(variantnames) * 3,
                 ncol = 4)
oddrows <- which(1:nrow(outtab) %% 3 == 1)
evenrows <- which(1:nrow(outtab) %% 3 == 2)
outtab[oddrows,-1] <- round(100 * decomptable)
outtab[evenrows,-1] <- inttabl
outtab[oddrows,1] <- variantnames
outtab <- as.data.frame(outtab)
names(outtab) <- c("Model", "Ideology", "Idiosyncrasy", "Instability")

outputtable <- data.frame('Model'=variantnames,
           Ideology= round(100*decomptable[,1]),
           Idiosyncrasy= round(100*decomptable[,2]),
           Instability= round(100*decomptable[,3]))
```

```{r, echo=FALSE, results = "asis"}
emphasize.italics.rows(evenrows)
pandoc.table(outtab, justify = c("right", rep("center", 3)))
```
